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Calculations in Quantum Chromodynamics are typically performed using a method 
pioneered by Ellis, Ross and Terrano in 1981. In this method, one combines nu- 
merical integrations over the momenta of final state particles with analytical inte- 
grations over the momenta of virtual particles. I discuss a more flexible method in 
which one performs all of the integrations numerically. 



1 Introduction 

This talk concerns a method for performing perturbative calculations in Quan- 
tum Chromodynamics and other quantum field theories, a The method may 
prove useful for cross sections in which one measures something about the 
hadronic final state, so that one cannot use the powerful but very special 
methods that apply to completely inclusive quantities like the total cross sec- 
tion for e + e~ — > hadrons. Examples of such cross sections include jet cross sec- 
tions in hadron-hadron and lepton-hadron scattering and in e + e~ — > hadrons. 

Many calculations of this kind have been carried out at next-to-leading 
order in perturbation theory, resulting in an impressive confrontation between 
theory and experiment. These calculations have been based on a method 
introduced by Ellis, Ross, and Terranoa in the context of e + e~ — > hadrons. 
Stated in the simplest terms, one performs some integrations over momenta pi 
analytically, others numerically. I shall argue that it is possible instead to do 
all of these integrations numerically. This has some advantages, principally in 
the flexibility that it allows. 

In this talk, I address specifically the calculation of three-jet-like infrared 
safe observables in e + e~ — > hadrons at next-to-leading order, that is order a 2 s . 
Examples of such observables include the thrust distribution, the fraction of 
events that have three jets, and the energy-energy correlation function. There 
are, of course, several existing computer programs that can calculate any such 
quantity, the earliest of these being that of Kunszt and Nason, El which was 
developed in 1989. What is new in this talk is not the result achieved, but the 
method of calculation. 



2 The quantity to be calculated. 



The order a 2 s contribution to the observable being calculated has the form 
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There are, of course, infrared divergences associated with Eq. (|l]); we under- 
stand that an infrared cutoff has been supplied. In Eq. (|lj), the daffl are the 
order a 2 s contributions to the parton level cross section, calculated with zero 
quark masses. Each contains momentum and energy conserving delta func- 
tions. The daffl include ultraviolet renormalization in the MS scheme. The 
functions S describe the measurable quantity to be calculated. Since_we cal- 
culate a "three-jet" quantity, 52 = 0. The measurement, as specified cl by the 
functions S n , is to be infrared safe: the S n are smooth functions of the parton 
momenta and 

S n+1 (kx, Xkn, (1 - X)k n ) = S n (ki, ...,k n ) (2) 

for < A < 1. That is, collinear splittings and soft particles do not affect the 
measurement. 

It is convenient to calculate a quantity that is dimensionless. Let us, 
therefore, define the functions S n so that they are dimensionless and eliminate 
the remaining dimensionality in the problem by dividing by do, the total e + e~ 
cross section at the Born level. We also remove the factor of (a s /ir) 2 . Thus, 
we calculate 

1= , , 3 



3 Setting up the calculation 

We note that I is a function of the cm. energy \fs and the MS renormalization 
scale li. We will choose ll to be proportional to -y/s: (i — Auv\fs- Then X 
depends on A. But, because it is dimensionless, it is independent of y/s. This 
allows us to write 
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(4) 



where h is any function with 



/ dyG h(yG) = 1. (5) 

Jo 

The integration over \fs eliminates the energy conserving delta function in X. 
This is important in allowing cancellations of infrared singularities between 
contributions with 3- and 4-particle intermediate states to take place. The 
physical meaning is that, by smearing in the energy y/s, we force the time 
variables in the two current operators that create the hadronic state to be 
within 1/s/s of each other. Thus we have a truly short distance problem. 



4 The Ellis, Ross, and Terrano method 

I now describe how one would calculate 1 using the Ellis- Ross- Terrano method. 
Each partonic cross section in Eq. (Q) can be expressed as an amplitude times 
a complex conjugate amplitude, as illustrated in Fig. [lj. One must calculate 
the amplitudes in 4 — 2e dimensions. (In the case of the process. e + e~ — * 
hadrons, this calculation was performed by Ellis, Ross, and Terrano. a) For tree 
diagrams, the calculation is straightforward. For loop diagrams, this involves 
an integration, which is performed analytically. The integrals are divergent 
in four dimensions, so one obtains divergent terms proportional to 1/e 2 and 
1/e in addition to terms that are finite as e — > 0. Having the amplitudes 
and complex conjugate amplitudes, one must now multiply by the functions 
S n and integrate over the final state parton momenta. These integrations are 
too complicated to perform analytically, so one must use numerical methods. 
Unfortunately, the integrals are divergent at e = 0. Thus one must split the 
integrals into two parts. One part can be divergent at e = but must be 
simple enough to calculate analytically. The other part can be complicated, 
but must be convergent at e = 0. One calculates the simple, divergent part 
and cancels the 1/e 2 and 1/e pole terms against the pole terms coming from 
the virtual loop diagrams. This leaves the complicated, convergent integration 
to be performed numerically. 

This method is a little bit cumbersome, but it works and has been enor- 
mously successful. However it has proven to be difficult to apply the method 
in the case of two virtual loops. Even with one virtual loop, the method is 
not very flexible. In modern implementions of the Ellis- Ross- Terrano method, 
one can easily choose the functions but any other modification of the in- 
tegrand requires one to recalculate the amplitudes, and the modification must 
be simple enough that one can calculate the amplitudes in closed form. 




Figure 1: Illustration of the method of Ellis, Ross and Terrano. The dots on the parton lines 
entering the final state represent the measurement function. Some of the contributions to 
the amplitudes are illustrated in the lower part of the figure. For four particle intermediate 
states, one multiplies tree amplitudes. For three particle intermediate states, one multiplies 
a tree amplitude by a loop amplitude. 

5 Calculation by numerical integration 

Let us, therefore, inquire whether there is any other way that one might per- 
form such a calculation. We note that the quantity X can be expressed in terms 
of cut Feynman diagrams, as in Fig. ^. Each diagram is a three loop diagram, 
so we have integrations over loop momenta ff^ and 



5.1 Performing the energy integrations 

We first perform the energy integrations. For the graphs in which four parton 
lines cross the cut, there are four mass-shell delta functions S(kj). These delta 
functions eliminate the three energy integrals over I®, ^2> an< ^ ^3 as weu as ^ ne 
integral (||) over \fs. For the graphs in which three parton lines cross the cut, 
we can eliminate the integration over y/s and two of the £°j integrals. One 
integral over the energy E in the virtual loop remains. The integrand contains 
a product of factors 

E-Q°j- luj +ie E-Q°j + uj-ie ' 1 ' 

where E — Q j is the energy carried by the Jth propagator around the loop 
and ujj is the absolute value of the momentum carried on that propagator. We 




Figure 2: Some cut Feynman diagrams that contribute to e+e — > hadrons. 

perform the integration by closing the integration contour in the lower half E 
plane. This leads to n terms for a virtual n point subgraph. In the Jth term, 
the propagator energy E — Qj is set equal to the corresponding uij and there 
is a factor l/(2u>j). Note that the entire process of performing the energy 
integrals amounts to some simple algebraic substitutions. 

5.2 Cancellations of pinch singularities 

Let us denote the contribution to T from the cut C of graph G by T{G, C). 
This contribution has the form 

AG, C) = J d 3 h d% d% g(G, C; £), (7) 

where I denotes the set of loop momenta {£1,(2, £3}- The functions g have some 
singularities, called pinch singularities, that cannot be avoided by deforming 
the integration contour and some non-pinch singularities that can be avoided 
by a contour deformation. 

Let us consider the pinch singularities, following the analysis of Sterman.B 
The pinch singularities occur when one parton branches into two partons with 
collincar momenta or when one parton momentum goes to zero. These sin- 
gularities can lead to logarithmic divergences in the corresponding integral. 
(There is a complication associated with the q^q v terms in gluon self-energy 
subgraphs but I ignore this complication in this talk.) If we were to calculate 



the total cross section by using measurement functions S n = 1 in g, then the 
singularities would cancel between the functions g associated with the various 
cuts C of the same graph G. The underlying reason is unitarity. Now in our 
case of measured shape variables, the values of S n corresponding to different 
cuts C are different. This would ruin the cancellation, except that just at the 
collinear or soft points the functions S n match. Thus the singularities present 
in the individual g(G,C;£) cancel in the sum, ^ c g(G, C; £). To be precise, 
the cancellation reduces the strength of the singularity from a strength suffi- 
cient to give a logarithmically divergent integral to one that gives a convergent 
integral. 

5.3 Avoiding the non-pinch singularities 

The function g also has singularities that can be avoided by deforming the 
integration contours, the non-pinch singularities. These singularities do not 
cancel, but the Feynman rules provide an ie prescription that tells us that 
we should deform the £ integration contour into the complex £ plane so as 
to avoid the singularity. Here deforming the contour means replacing £ by a, 
complex vector £ + in. Then one simply chooses the imaginary part, k, of the 
loop momentum as a function of the real part, £, and supplies the appropriate 
jacobian J. Then 



1(G, C)= d% d 3 £ 2 d 3 £ 3 J{G, C; £) g{G, C;£ + in(G, C; £)). (8) 



1(G) = / d 3 h d 3 £ 2 d 3 £ 3 J( G > C '> l ) f( G > C '> £ + iK ( G > C ' £ )) ( 9 ) 



is finite. One can simply compute it by Monte Carlo integration after removing 
from the integration tiny regions near the collinear and soft singular points, 
where roundoff errors spoil the cancellation of the individual contributions. 

Note the significance of putting the summation over cuts inside the inte- 
gral. When we sum over cuts for a given point in the space of loop momenta, 
the soft and collinear divergences cancel because the cancellation is built into 
the Feynman rules. If we were to put the sum over cuts outside the integra- 
tion, as in the Ellis-Ross- Terrano method, then the individual integrals would 
be divergent. The calculation would thereby be rendered more difficult. 




5.4 The integration 



With the contours appropriately chosen, the integral 



C 



Figure 3: Two cuts of a diagram with a self-energy subgraph. 



6 Additional points 

Eq. (||) represents the main point of this talk. There are two additional points 
of principle that I discuss below: renormalization and the special treatment 
required for self-energy subdiagrams. 



6.1 Renormalization 

Some of the virtual loop subgraphs are ultraviolet divergent. Normally, one 
would renormalize these subgraphs by performing loop integrations in 4 — 2e 
space-time dimensions and subtracting the resulting pole term, c/e. Clearly 
that is not appropriate in a numerical integration. However, one can subtract 
instead an integral in 4 space-time dimensions such that, in the region of large 
loop momenta, the integrand of the subtraction term matches the integrand of 
the subdiagram in question. The integrand of the subtraction term can depend 
on a mass parameter fj, in such a way that the subtraction term has no infrared 
singularities. (For example, in the simplest case one can have denominators of 
the form + ft 2 ].) Then, one can easily arrange the definition so that this 
subtraction has exactly the same effect as the conventional MS subtraction 
with scale parameter /x. 



6.2 Self- energy subgraphs 

Virtual self-energy subgraphs, as in the right hand diagram in Fig. [| require 
a special treatment. Consider a quark self-energy subdiagram — zS with one 
adjoining virtual propagator and one adjoining cut propagator. This combina- 
tion really represents a field strength renormalization for the quark field, and 
is interpreted as 

1 



(10) 



In order to take the q 2 — ► limit here while at the same time maintaining the 
cancellation of collinear divergences, we write 

m<i)i _ a 2 c f d3i g \k + h -(\k + \ + \k-\)k + -j 
i 2 ( 2 ^) 3 F J \k + \\U[{\k + \ + \U) 2 -{q G ) 2 ^ 

where k± = \q ± I. This expression is obtained by writing the left-hand 
side as a dispersive integral with the cut self-energy graph appearing as the 
discontinuity. When the virtual self-energy is written in this form, the q 2 — > 
limit is smooth and, in addition, the cancellation between the two cut diagrams 
in Fig. U in the collinear limit £ oc q is manifest. It should be noted that the 
integral in Eq. ( |TT| ) is ultraviolet divergent and requires renormalization, which 
can be performed with an ad hoc subtraction as described above. 

One may expect that the representation of virtual propagator corrections 
in terms of the cut propagator will prove to be convenient in future modifica- 
tions of the method described here. One may want to make modifications to 
the gluon propagator, in particular, in order to implement a running strong 
coupling and to insert models for the long distance propagation of gluons. In 
addition, one may want to modify propagators so that partons can enter the 
final state with q 2 > so as to mesh an order a 2 perturbative calculation 
with a parton shower Monte Carlo program. In either case, the modified vir- 
tual gluon propagator should be related by a dispersive representation to the 
modified final state created by an outgoing gluon. 



7 An example at order a s 

I have argued that a completely numerical calculation of the e + e~ observables 
considered here is possible in principle. The question now arises whether the 
principle can be put into operation in an actual computer program. 

Consider first a simple example at one lower order in a s : the total cross 
section for e + e" — > hadrons at order a s . The cut diagrams are shown in 
Fig. |. 

Let us define to be the sum of the top graphs in Fig. |] and ^(/i) 

to be the sum of the bottom graphs, so that 

* = a Q {l + {A 1 (ri+A 2 (ti)}^}. (12) 

Then a simple analytical calculation shows that 



Ai(/i)=-|-|ln( At />/i), A 2 ( M ) = ~ + ^m( / VVi)- (13) 




Figure 4: The cut diagrams for the total cross section for e + e — ► hadrons at order a s 

I have constructed a demonstration computer program to calculate A\ (/i) 
and Ai{p) numerically along the lines described above. This example has 
UV renormalization, propagators to be treated dispersively, singularities to 
be avoided by deforming the contour, soft singularities that must cancel, and 
collinear singularities that must cancel. Thus it is not trivial from a numerical 
point of view. The results are shown in Fig. |^. Evidently, the method works. 



8 An example at order a 2 

In order to further test the methods described here, I have constructed a 
demonstration computer program El to calculate three-jet-like observables for 
e + e~ annihilation at order a 2 . I have used the program to calculate da^ /dT, 
the order a 2 contribution to the thrust distribution. More precisely, I have cal- 
culated the ratio R(T) of d<M /dT to [da^ /dT] KN , where where [da® /d^W 
is a fit to the tabulated results for da^ /dT as given by Kunszt and Nason.Q In 
the range 0.71 < T < 0.95, the function [da^ /dT) KN varies by about a factor 
of 30. The ratio R(T) should be 1. The results are reported in Fig. [| Again, 
it appears that the numerical method of calculation works in a practical sense. 



9 Outlook 

I am working to test the code of the demonstration program. Indeed, between 
the time of the RADCOR98 conference and the time of this writing, I have 
found and corrected a small error in the program logic which resulted in an 
error in the calculation of about 1%. With this bug corrected, the results 
produced by the current version of the program are better than those shown 
in Fig. ||. After testing the code, I will have to document the code and publish 
a detailed description of the algorithms used. 

With this code in hand, or with improved code from other authors, I 
anticipate attacks on more difficult problems than the problem discussed here. 




Figure 5: Result of the numerical calculation of the two pieces of the total cross section for 
e + e — — > hadrons at order a s . The numerically calculated functions and j4.2(m) are 

plotted versus the renormalization scale. The lines are the analytical results. 
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Figure 6: Ratio R of the order contribution to the thpjist distribution as calculated here 
to the same quantity as calculated by Kunszt and Nason.U The horizontal lines represent the 
expected result, 1, with an error band based on the quoted errors in the results of Kunszt 
and Nason and the uncertainties in fitting a smooth function to these results. The error 
bars on the computed points represent statistical errors. There is some correlation expected 
between neighboring points. 



It remains to be seen for what problems the completely numerical method 
will prove to be more powerful than the analytical/numerical method that has 
served us so well up to now. 
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